Show code
library(dplyr)
library(readr)
library(ggplot2)
library(tigris)
library(sf)
library(DT)Caitlin Uang
December 15, 2025
To acquire the national fast food count per state and county, this project will use the NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990-2021 dataset and filter year 2019 for more accurate pre-COVID numbers. The dataset at census tract level is used.
# Cache shapefiles locally (faster repeat runs)
options(tigris_use_cache = TRUE)
# load dataset
eatdrink_tract <- read_csv("nanda_eatdrink_Tract20_1990-2021_01P.csv")
fastfood_raw <- eatdrink_tract |>
filter(year == 2019) |>
arrange(desc(count_fastfood)) |>
select(tract_fips20,totpop,count_fastfood,den_fastfood,)
tibble(fastfood_raw)# A tibble: 85,528 × 4
tract_fips20 totpop count_fastfood den_fastfood
<chr> <dbl> <dbl> <dbl>
1 06073008511 4989 22 4.41
2 12086009010 6716 22 3.28
3 06059062610 10440 20 1.92
4 12023110500 7764 20 2.58
5 17031839100 7545 18 2.39
6 04013216816 6656 17 2.55
7 12095017001 4652 16 3.44
8 06071002109 6121 15 2.45
9 29183312400 5089 15 2.95
10 40109108604 3664 15 4.09
# ℹ 85,518 more rows
pop_2019 <- read_csv("ACSDT5Y2019.B01003-Data.csv")
pop_clean <- pop_2019 |>
mutate(tract_fips20 = substr(GEO_ID, nchar(GEO_ID) - 10, nchar(GEO_ID)))
fastfood_tract <- fastfood_raw |>
left_join(pop_clean, by = "tract_fips20") |>
select(tract_fips20,GEO_ID,totpop,B01003_001E,den_fastfood,count_fastfood) |>
rename(acs_pop = B01003_001E) |>
mutate(acs_pop = as.numeric(acs_pop)) |>
filter(acs_pop > 0)
datatable(fastfood_tract)#State look up
state_lookup <- states(year = 2019) |>
st_drop_geometry() |>
select(STATEFP, STUSPS, NAME) |>
rename(
State_Abbr = STUSPS,
State_Name = NAME
)
# Aggregate tract to county
den_fastfood_county <- fastfood_tract |>
mutate(county_fips20 = substr(tract_fips20, 1, 5)) |>
group_by(county_fips20) |>
summarize(
total_fastfood = sum(count_fastfood, na.rm = TRUE),
total_pop = sum(acs_pop, na.rm = TRUE),
den_fastfood = (total_fastfood / total_pop) * 1000,
.groups = "drop"
) |>
filter(total_pop > 0)
# Load shapefile and join
county_shapes <- counties(year = 2019, cb = TRUE) |>
mutate(GEOID = as.character(GEOID))
fastfood_map_county <- county_shapes |>
left_join(den_fastfood_county, by = c("GEOID" = "county_fips20")) |>
left_join(state_lookup, by = "STATEFP") |>
mutate(
den_fastfood_plot = ifelse(den_fastfood == 0, NA, den_fastfood)
)
# Create datatable
datatable(
fastfood_map_county |>
st_drop_geometry() |>
arrange(desc(den_fastfood)) |>
select(
County = NAME,
State = State_Abbr,
State_Name,
total_fastfood,
total_pop,
den_fastfood
),
caption = "Fast Food Restaurants per 1,000 People (County, 2019)",
options = list(pageLength = 10, scrollX = TRUE)
) |>
formatRound("total_fastfood", 0, mark = ",") |>
formatRound("den_fastfood", 2)fastfood_map_county |>
st_drop_geometry() |>
ggplot(aes(x = den_fastfood)) +
geom_histogram(
bins = 30,
fill = "#9badff",
color = "white",
linewidth = 0.3
) +
coord_cartesian(xlim = c(0, 1)) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "Distribution of Fast Food Density per County, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)",
x = "Fast Food Density",
y = "# of Counties"
) +
theme(
axis.line = element_line(color = "black", linewidth = 0.4),
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0),
panel.background = element_rect(fill = "white")
)ggplot(fastfood_map_county) +
geom_sf(aes(fill = den_fastfood), color = NA, linewidth = 0.1) +
scale_fill_gradientn(
colours = c("grey90", "#e792e7", "#5b7cff", "#0432FF"),
name = "Fast Food\nDensity",
limits = c(NA, NA),
na.value = "#ffffff",
trans = "sqrt"
) +
coord_sf(
xlim = c(-125, -66),
ylim = c(24, 50),
expand = FALSE
) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "County, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0)
)# Aggregate tract to state
fastfood_state <- fastfood_tract |>
mutate(state_fips20 = substr(tract_fips20, 1, 2)) |>
group_by(state_fips20) |>
summarize(
total_fastfood = sum(count_fastfood, na.rm = TRUE),
total_pop = sum(acs_pop, na.rm = TRUE),
den_fastfood = (total_fastfood / total_pop) * 1000,
.groups = "drop"
) |>
filter(total_pop > 0)
# Load state shapefile and join
state_shapes <- states(year = 2019, cb = TRUE) |>
mutate(STATEFP = as.character(STATEFP))
fastfood_map_state <- state_shapes |>
left_join(fastfood_state, by = c("STATEFP" = "state_fips20"))
# Creat datatable
datatable(
fastfood_map_state |>
st_drop_geometry() |>
arrange(desc(den_fastfood)) |>
slice_max(den_fastfood, n = 10) |>
select(State = NAME, den_fastfood,total_pop,total_fastfood),
rownames = FALSE,
caption = "Fast Food Restaurants per 1,000 People (State, 2019)",
options = list(
searching = FALSE,
scrollX = FALSE,
autoWidth = TRUE
)
) |>
formatRound("den_fastfood", 2)# Create heat map
ggplot(fastfood_map_state) +
geom_sf(aes(fill = den_fastfood), color = "white") +
scale_fill_gradient(
low = "#ffffff",
high = "#932092",
name = "Fast Food\nDensity", limits = c(NA,NA)
)+
coord_sf(
xlim = c(-125, -66),
ylim = c(24, 50),
expand = FALSE
) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "State, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0)
)# Rank and flag top 10 states
fastfood_map_state_topflag <- fastfood_map_state |>
arrange(desc(den_fastfood)) |>
mutate(
rank = row_number(),
den_top10 = ifelse(rank <= 10, den_fastfood, NA_real_)
)
ggplot(fastfood_map_state_topflag) +
geom_sf(aes(fill = den_top10),
color = "white",
linewidth = 0.3) +
scale_fill_distiller(
palette = "Blues", # Nice clean blues
direction = 1, # Light → dark
na.value = "#e0e0e0", # Light grey for non–top 10
name = "Fast-Food per\n1,000 People"
) +
labs(
title = "Top 10 States with the Most Fast-Food Restaurants per 1,000 People",
subtitle = "Non–Top 10 States Shown in Grey",
caption = "Source: NaNDA Tract20 (1990–2021) & U.S. Census ACS 2019"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold"),
plot.subtitle = element_text(size = 12),
plot.caption = element_text(size = 9, hjust = 0),
plot.margin = margin(10, 10, 10, 10)
)